R Script

## Loading and installing packages

knitr::opts_chunk$set(echo = F, message = F, warning = F)

packages <- c("rgdal", "foreign", "gdata", "ggmap", "ggplot2",
              "plyr", "rgeos", "sf", "ggrepel", "dplyr", "sp", "aspace",
              "spdep", "bookdown", "stringr", "maptools", "leaflet", "broom", "stargazer",
              "RColorBrewer", "mctest", "olsrr", "knitr", "spgwr", "htmlwidgets", "htmltools")

package.check <- lapply(packages, FUN = function(x) {
  if (!require(x, character.only = T)) install.packages(x)
  if (! (x %in% (.packages() )))  library(x, character.only = T)
})

# leaflet map titles
tag.map.title <- tags$style(HTML("
  .leaflet-control.map-title { 
    transform: translate(-50%,20%);
    position: fixed !important;
    left: 50%;
    text-align: center;
    padding-left: 10px; 
    padding-right: 10px; 
    background: rgba(255,255,255,0.75);
    font-weight: bold;
    font-size: 28px;
  }
"))



p <- read.csv("../Research Data/2010_Census_Populations_by_Zip_Code.csv")
load("../Research Data/crime.RData")
load(file = "../Research Data/zipCrimes.RData")

names(zc)[2] <- names(p)[1] <- "zip"

zCrimes <- zc
zCrimes@data <- merge(zc@data, p, by = "zip")
# zCrimes <- zCrimes[!is.na(zCrimes$crimes),]

# writeOGR(obj=zCrimes, driver="ESRI Shapefile", "../Research Data/zCrimes")



par(mfrow=c(1,3))
hist(zCrimes$crimes, xlab = "crimes since 2010", main=NULL)
hist(zCrimes$Average.Household.Size, xlab = "average household size in 2010", main=NULL)
hist(zCrimes$Median.Age, xlab = "median age in 2010", main=NULL)
par(mfrow=c(1,1))


par(mfrow=c(1,2))
plot(zCrimes$Median.Age, log(zCrimes$crimes), xlab = "median age", ylab = "ln(crimes)")
plot(zCrimes$Average.Household.Size, log(zCrimes$crimes), , xlab = "average household size", ylab = "ln(crimes)")
par(mfrow=c(1,1))


centroids <- as.data.frame(gCentroid(zCrimes,byid=TRUE))

# pal <- colorNumeric(
#   palette = "YlGnBu",
#   domain = zCrimes$crimes
# )

qpal <- colorQuantile("YlGnBu", zCrimes$crimes, n = 5)
qHH <- colorQuantile("YlOrRd", zCrimes$Average.Household.Size, n = 5)

leaflet(zCrimes) %>% addPolygons(weight = 1, smoothFactor = 0.5,
    opacity = 1.0, fillOpacity = 0.5,
    color = ~qpal(crimes),
    highlightOptions = highlightOptions(color = "white", weight = 2,
      bringToFront = TRUE)) %>% 
                addCircles(lng = ~centroids$x, lat = ~centroids$y, weight = 1, color = ~qHH(Average.Household.Size),
                  radius = 800, popup = ~zCrimes$zip, opacity = 0.9, fillOpacity = 0.8 ) %>% 
                addTiles() %>%
    addLegend(pal = qpal, values = ~crimes, opacity = 1) %>%
    addLegend(pal = qHH, values = ~Average.Household.Size, opacity = 1)

qMA <- colorQuantile("YlOrRd", zCrimes$Median.Age, n = 5)


leaflet(zCrimes) %>% addPolygons(weight = 1, smoothFactor = 0.5,
    opacity = 1.0, fillOpacity = 0.5,
    color = ~qpal(crimes),
    highlightOptions = highlightOptions(color = "white", weight = 2,
      bringToFront = TRUE)) %>% 
                addCircles(lng = ~centroids$x, lat = ~centroids$y, weight = 1, color = ~qMA(Median.Age),
                  radius = 800, popup = ~zCrimes$zip, opacity = 0.9, fillOpacity = 0.8 ) %>% 
                addTiles() %>%
    addLegend(pal = qpal, values = ~crimes, opacity = 1) %>%
    addLegend(pal = qMA, values = ~Median.Age, opacity = 1)


zCrimes@data$id = rownames(zCrimes@data)
crimePoints = fortify(zCrimes, region="id")
crimesDf = join(crimePoints, zCrimes@data, by="id")


# c$weapon <- !c$Weapon.Description %in% c("", "STRONG-ARM (HANDS, FIST, FEET OR BODILY FORCE)", "VERBAL THREAT")
# c$verbal <- c$Weapon.Description %in% "VERBAL THREAT"
# c$fists <- c$Weapon.Description %in% "STRONG-ARM (HANDS, FIST, FEET OR BODILY FORCE)"
# 
# clean_coords <- gsub(pattern = '[()]', replacement = '', x = c$Location)
# split_coords <- str_split(string = clean_coords, pattern = ', ', n = 2, simplify = T)
# c$lat <- as.numeric(split_coords[,1])
# c$lon <- as.numeric(split_coords[,2])


zCrimesLoc <- cbind(zCrimes@data, centroids)

ggplot(crimesDf) + 
  geom_polygon(aes(long,lat,group=group, fill = crimes)) +
  coord_equal() + 
  stat_ellipse(data = subset(zCrimesLoc, Median.Age > 40), aes(x = x, y = y), level=0.5, color = "red") +
  geom_point(data = subset(zCrimesLoc, Median.Age > 40), aes(x = mean(x), y = mean(y)), color = "red", size = 0.5) +
  stat_ellipse(data = subset(zCrimesLoc, Median.Age <= 40), aes(x = x, y = y), level=0.5, color = "orange") +
  geom_point(data = subset(zCrimesLoc, Median.Age <= 40), aes(x = mean(x), y = mean(y)), color = "orange", size = 0.5) +
  theme_void()


ggplot(crimesDf) + 
  geom_polygon(aes(long,lat,group=group, fill = crimes)) +
  coord_equal() + 
  stat_ellipse(data = subset(zCrimesLoc, Average.Household.Size > 3), aes(x = x, y = y), level=0.5, color = "red") +
  geom_point(data = subset(zCrimesLoc, Average.Household.Size > 3), aes(x = mean(x), y = mean(y)), color = "red", size = 0.5) +
  stat_ellipse(data = subset(zCrimesLoc, Average.Household.Size <= 3), aes(x = x, y = y), level=0.5, color = "orange") +
  geom_point(data = subset(zCrimesLoc, Average.Household.Size <= 3), aes(x = mean(x), y = mean(y)), color = "orange", size = 0.5) +
  theme_void()


regHH <- lm(log(crimes) ~ Average.Household.Size + Median.Age, data = zCrimes)
stargazer(regHH, no.space = T, type = "html")

# Multicollinearity

mc <- omcdiag(zCrimes@data[colnames(zCrimes@data) %in% c("Average.Household.Size", "Median.Age")], log(zCrimes$crimes))

kable(mc$odiags, "html", booktabs = T)


# HETEROSKEDASTICITY

# ols_test_bartlett(as.data.frame(zCrimes@data$"Average.Household.Size"), as.data.frame(zCrimes@data$"Median.Age"))
bp <- ols_test_breusch_pagan(regHH)

m <- matrix(c(bp$bp, bp$p, 1))
rownames(m) <- c("Chi2", "Prob > Chi2", "DF")
kable(m, "html", booktabs = T, row.names = T, caption = "Breusch-Pagan test")


# Spatial Dependence

spatmatrix <- poly2nb(zCrimes)
w <- nb2listw(spatmatrix, zero.policy = T)

# Moran's I
moran <- lm.morantest(regHH, listw = w, zero.policy = T)# Moran’s I test for residual spatial autocorrelation
moranStat <- matrix(c(moran$statistic, moran$p.value, 1))
rownames(moranStat) <- c("statistic", "p.value", "parameter")
colnames(moranStat) <- "Moran's I (error)"

# Lagrange Multiplier test
LM <- lm.LMtests(regHH, listw = w, zero.policy = T, test = c("LMerr", "RLMerr", "LMlag", "RLMlag", "SARMA")) 

LMrows <- c("statistic", "p.value", "parameter")
m <- matrix(NA, nrow = 3, ncol = length(LM))
rownames(m) <- LMrows
colnames(m) <- c("Lagrange Multiplier (error)", "Robust LM (error)", "Lagrange Multiplier (lag)", "Robust LM (lag)", "Lagrange Multiplier (SARMA)")
# colnames(m) <- 
  
# everything in one table
for(i in 1:length(LM)){
  r <- 1
  for(j in LMrows){
    m[r, i] <- LM[[i]][[j]]
    r <- r + 1
  }
}
m <- cbind(moranStat, m)
kable(m, "html", booktabs = T, row.names = T)



# Residuals

#extract residuals
zCrimes@data$lmRes[!is.na(zCrimes$crimes)] <- resid(regHH) #residual lm

#view residuals for linear model
qpal <- colorQuantile("OrRd", zCrimes@data$lmRes, n=9) 

leaflet(zCrimes) %>%
  addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2, color = ~qpal(lmRes)
  ) %>%
  addTiles()


#run a spatial lag model
lagModel <- lagsarlm(log(crimes) ~ Average.Household.Size + Median.Age, data = zCrimes, w,
                  zero.policy = T)
stargazer(lagModel, no.space = T, type = "html")
#note that rho is listed below other model coefs.


#run a spatial error model
errModel <- errorsarlm(log(crimes) ~ Average.Household.Size + Median.Age, data = zCrimes, w, zero.policy = T)
stargazer(errModel, no.space = T, type = "html")




#extract residuals for lag model
zCrimes@data$lagRes[!is.na(zCrimes$crimes)] <- resid(lagModel) #residuals lag

#view residuals for lag model
qpal<-colorQuantile("OrRd", zCrimes@data$lagRes, n=9) 

title <- tags$div(tag.map.title, HTML("spatial lag model residuals"))  

leaflet(zCrimes) %>%
  addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2, color = ~qpal(lagRes)
  ) %>%
  addTiles() %>%
  addControl(title, position = "topleft", className="map-title")


#extract residuals for err model
zCrimes@data$errRes[!is.na(zCrimes$crimes)] <- resid(errModel) #residual err

#compare with residuals for err model
qpal<-colorQuantile("OrRd", zCrimes@data$errRes, n=9) 

title <- tags$div(tag.map.title, HTML("spatial error model residuals"))  

leaflet(zCrimes) %>%
  addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2, color = ~qpal(errRes)
  ) %>%
  addTiles() %>%
  addControl(title, position = "topleft", className="map-title")



centr <- gCentroid(zCrimes, byid = TRUE)  #byid tells function to add centroids per block group

# create SpatialPointsDataFrame to add all data from chi.poly@data to centr@data
centr <- SpatialPointsDataFrame(centr, data= zCrimes@data) 


#Begin GWR prep: First, calculate bandwidth

#calculate optimal bandwidth for model (uses cross validation to find kernel bandwidth that
#generates best model...e.g.-kernel that minimizes RSS)

GWRbandwidth <- gwr.sel(log(crimes) ~ Average.Household.Size + Median.Age, data=centr, verbose = F, 
                        coords=coordinates(centr),adapt=T)# coordinates() extracts x, y coordinates
#adapt tells function to calculate adaptive bandwidth

#run the gwr model
gwr.model <- gwr(log(crimes) ~ Average.Household.Size + Median.Age, 
                 data=subset(centr, !is.na(crimes)), coords=coordinates, adapt=GWRbandwidth,
                 hatmatrix=TRUE, se.fit=TRUE)
#print the results of the model

#extract results for each variable
results<-as.data.frame(gwr.model$SDF)

#attach coefficients to original dataframe
centr$Median.AgeCoef[!is.na(zCrimes$crimes)]  <- results$Median.Age
centr$Average.Household.SizeCoef[!is.na(zCrimes$crimes)]  <- results$Average.Household.Size

#now plot the various GWR coefficients

#Does the median age influence violent crime uniformly across LA?

title <- tags$div(tag.map.title, HTML("median age influence on violent crime"))  

qpal<-colorQuantile("OrRd", centr@data$Median.Age, n=9) 

leaflet(centr) %>%
  addCircleMarkers(radius=3,color = ~qpal(centr@data$Median.Age)
  ) %>%
  addTiles() %>%
  addControl(title, position = "topleft", className="map-title")



#Does the average household size influence violent crime uniformly across LA?
title <- tags$div(tag.map.title, HTML("average household size influence on violent crime"))  

qpal<-colorQuantile("OrRd", centr@data$Average.Household.Size, n=9) 

leaflet(centr) %>%
  addCircleMarkers(radius=3,color = ~qpal(centr@data$Average.Household.Size)
  ) %>%
  addTiles() %>%
  addControl(title, position = "topleft", className="map-title")